function k = smoothkern(f, t)
% smoothkern  - build N-dim smoothing kernel from FWHM
%
% FORMAT:       k = smoothkern(f [, t])
%
% Input fields:
%
%       f           1xF FWHM values (e.g. [2, 2, 2])
%       t           optional weight threshold (default: eps ^ (1/numel(f)))
%
% Output fields:
%
%       k           kernel for convolution

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
% argument check
if nargin < 1 || ...
   ~isa(f, 'double') || ...
    isempty(f) || ...
    numel(f) ~= max(size(f)) || ...
    ndims(f) > 2 || ...
    any(isinf(f) | isnan(f) | f < 0 | f > 32)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid argument.' ...
    );
end
f = f(:)';
if nargin < 2 || ...
   ~isa(t, 'double') || ...
    numel(t) ~= 1 || ...
    isinf(t) || ...
    isnan(t) || ...
    t < 0 || ...
    t > ((1/4) ^ numel(f))
    t = eps ^ (1 / numel(f));
end

% lower bound on FWHM is 0.29, for which the kernel function has no 
f = max(f, 0.29 * ones(size(f)));

% FWHM -> Gaussian parameter
f = f / sqrt(8 * log(2));

% max dist and dimension
md = round(6 * f);
fd = 2 * md + 1;

% kernel
k = ones([fd, 1]);

% for each dim multiply with linear kernel
for nd = 1:numel(f)
    xd = fd;
    rd = fd;
    xd(nd) = 1;
    rd([1:nd-1, nd+1:end]) = 1;
    ed = exp(- (-md(nd):md(nd)) .^ 2 ./ (2 * f(nd) .^ 2));
    ed = ed ./ sum(ed);
    k = k .* repmat(reshape(ed, [rd, 1]), [xd, 1]);
end

% threshold
if t > 0
    k(k < t) = 0;
    
    % give back smallest possible kernel
    s = cell(1, numel(f));
    for sc = 1:numel(s)
        ks = k;
        for s2 = 1:numel(s)
            if s2 ~= sc
                ks = sum(ks, s2);
            end
        end
        s{sc} = find(ks(:) > 0);
    end
    k = subsref(k, struct('type', '()', 'subs', {s}));
end

% make sure all elements sum to 1
k = k ./ sum(k(:));
